﻿#include "fft.h"

complex x[1024], *W;   //定义输入序列和旋转因子
double x_abs[1024];
int data_size = 0;   //定义数据长度
double PI = 4.0 * atan(1); //定义π 因为tan(π/4)=1 所以arctan（1）*4=π，增加π的精度

/**************************************复数******************************************/
// 复数运算
void add(complex a, complex b, complex* c)//定义结构体a、b和指针c    加法
{
    c->real = a.real + b.real; //c的目的是取出结构体中的数据
    c->imag = a.imag + b.imag;
}
void sub(complex a, complex b, complex* c)   //减法
{
    c->real = a.real - b.real;
    c->imag = a.imag - b.imag;
}
void mul(complex a, complex b, complex* c)   //乘法
{
    c->real = a.real * b.real - a.imag * b.imag;
    c->imag = a.real * b.imag + a.imag * b.real;
}

// 复数输出
void output(void)
{
    int i;
    for (i = 0; i < data_size; i++)
    {
        qDebug("%lf", x[i].real); //输出复数的实部
        if (x[i].imag >= 0.0001)
        {
            qDebug("+%lfj\n", x[i].imag);  //当复数的虚部大于0.0001时，输出+ 虚部 j的形式
        }
        else if (fabs(x[i].imag) < 0.0001)
        {
            qDebug("\n");//当虚部小于0.001时，跳过虚部，不输出
        }
        else
        {
            qDebug("%lfj\n", x[i].imag);//上述两个条件除外的形式，输出 虚部 j的形式
        }
    }
}


/*****************************************旋转因子*******************************************/
void transform(void)
{
    int i;
    W = (complex*)malloc(sizeof(complex) * data_size);
    for (i = 0; i < data_size; i++)
    {
        W[i].real = cos(2 * PI / data_size * i);
        W[i].imag = -1 * sin(2 * PI / data_size * i);
    }
}

/*******************************************算法实现******************************************/
// 数据点经过log(N) / log2级分割后重新排序
void change(void)
{
    complex temp;
    unsigned short i = 0, j = 0, k = 0;
    double t;
    for (i = 0; i < data_size; i++)
    {
        k = i;
        j = 0;
        t = (log(data_size) / log(2));  //算出序列的级数
        while ((t--) > 0)  //利用按位与以及循环实现码位颠倒
        {
            j = j << 1;
            j |= (k & 1);
            k = k >> 1;
        }
        if (j > i)    //将x(n)的码位互换
        {
            temp = x[i];
            x[i] = x[j];
            x[j] = temp;
        }
    }
    //output();
}

// 蝶形运算实现
void fft(void)
{
    int i = 0, j = 0, k = 0, m = 0;
    complex q, y, z;
    change();
    for (i = 0; i < log(data_size) / log(2); i++)  //蝶形运算的级数
    {
        m = 1 << i;   //移位 每次都是2的指数的形式增加，其实也可以用m=2^i代替
        for (j = 0; j < data_size; j += 2 * m)  //一组蝶形运算，每一组的蝶形因子乘数不同
        {
            for (k = 0; k < m; k++)  //蝶形结点的距离  一个蝶形运算 每个组内的蝶形运算
            {
                mul(x[k + j + m], W[data_size * k / 2 / m], &q);
                add(x[j + k], q, &y);
                sub(x[j + k], q, &z);
                x[j + k] = y;
                x[j + k + m] = z;
            }
        }
    }
}


void TXT_Read(void)
{
    QFile file(import_path);
    if(file.open(QIODevice::ReadOnly))
    {
        QTextStream stream(&file);
        while(!file.atEnd())
        {
            float buf;
            QStringList list=stream.readAll().split("  ");
            QListIterator<QString> li(list);
            for(int i=0;li.hasNext();i++)
            {
                buf=li.next().toDouble();
                x[i].real = buf;
            }
        }
    }
}

void TXT_Write(void)
{
    QFile file(output_path+"FFTOutput.txt");
    QFile::remove(output_path+"FFTOutput.txt");
    if(! file.open(QIODevice::Append|QIODevice::Text))  //append追加，不会覆盖之前的文件
    {
        qDebug("无法写入\r\n");
    }
    QTextStream out(&file);//写入
    for(int i=0;i<data_size;i++)
    {
        Calc_abs(i);
        out << x_abs[i]<< "  ";
    }
    file.close();
}
// 计算幅度
void Calc_abs(int i)
{
   x_abs[i] = sqrt(x[i].real * x[i].real + x[i].imag * x[i].imag);
}


//int main()
//{
//TXT_Read();
//transform();//变换序列顺序
//fft();//蝶形运算
//printf("输出FFT后的结果\n");
//output();//输出结果
//return 0;
//}
